Regularization of spherical and axisymmetric evolution codes in numerical relativity 
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Several interesting astrophysical phenomena are symmetric with respect to the rotation axis, 
like the head-on collision of compact bodies, the collapse and/or accretion of fields with a large 
variety of geometries, or some forms of gravitational waves. Most current numerical relativity 
codes, however, cannot take advantage of these symmetries due to the fact that singularities in the 
adapted coordinates, either at the origin or at the axis of symmetry, rapidly cause the simulation to 
crash. Because of this regularity problem it has become common practice to use full-blown Cartesian 
three-dimensional codes to simulate axi-symmetric systems. In this work we follow a recent idea 
of Rinne and Stewart and present a simple procedure to regularize the equations both in spherical 
and axi-symmetric spaces. We explicitly show the regularity of the evolution equations, describe 
the corresponding numerical code, and present several examples clearly showing the regularity of 
our evolutions. 

PACS numbers: 04.20.Ex, 04.25. Dm, 95.30.Sf, 



I. INTRODUCTION 

After forty years of research, the black hole collision 
problem can finally be considered solved. Though there 
are certainly still many details to be worked out, the re- 
sults from Pretorius ^] , the Brownsville and the Goddard 
groups 0, HI, and other groups that have followed, show 
that it is now possible to follow the numerical evolution 
of two black holes for several orbits, through the merger 
and subsequent ringing of the final merged black hole. 

Such tremendous progress in full three-dimensional 
(3D) numerical relativity, however, does not imply that 
there are no more interesting astrophysical situations 
that can be studied in spherical or axial symmetry, such 
as for example gravitational collapse or accretion onto 
compact objects. Accurate 3D simulations still require 
large computational resources, so that exploiting the ex- 
isting symmetries should allow important savings in com- 
putational time. Using coordinates adapted to the sym- 
metry, the number and complexity of the evolution equa- 
tions are reduced and thus the computational cost is also 
reduced. 

Nevertheless, the development of general purpose 
spherical/axi-symmetric codes in numerical relativity has 
been hampered by the lack of a generic method to deal 
with the singularities associated with the symmetry- 
adapted coordinate systems. For example, in the spheri- 
cally symmetric case described with spherical coordinates 
{r,9,(j)), the coordinates become singular at the origin 
r — Q. This implies that several terms in the evolution 
equations diverge as l/r, and even though local flatness 
guarantees that analytically all those terms should can- 
cel, such exact cancellation usually fails to hold in the 
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numerical description. A similar problem arises in sys- 
tems with axial symmetry when approaching the axis of 
symmetry. 

Several methods to deal with this problem have been 
proposed in the past. For example, one can choose spe- 
cific gauges that either eliminate or ameliorate the reg- 
ularity problem such as the areal (or radial) gauge in 
spherical symmetry, where the radial coordinate r is cho- 
sen in such a way that the proper area of spheres of con- 
stant r is always Anr'^. Similarly, in axial symmetry one 
can use the shift vector to guarantee that some metric 
components always vanish thus reducing the problem of 
regularity at the axis (for details see e.g. dlQ). Fur- 
thermore, there has been a lot of work on the construc- 
tion of axial codes that ensure that the metric remains 
smooth on the axis. For example, Garfinkle and Duncan 
describe in 7j a method that consists on the introduction 
of auxiliary variables which allow one to impose all the 
required regularity conditions on the extrinsic curvature. 
However, this method requires to solve, on every time 
slice, an elliptic equation for the lapse, the shift com- 
ponents and the conformal factor. A similar algorithm 
was presented by Choptuik et al. in [8], but adapted 
to the (2 -I- 1) -I- 1 formulation. Recently, another reg- 
ularization procedure was described in detail for the Z4 
system [ol, , by Rinne and Stewart in [ll|, [l^l , again 
also adapted to the (2 -I- 1) -I- 1 formulation. A different 
idea is the so-called "Cartoon method" , which consists in 
evolving three adjacent planes in Cartesian coordinates 
and then performing a tensor rotation to obtain bound- 
ary conditions p^ . However, as this method uses a tridi- 
mensional code, it is still more computationally more ex- 
pensive than an axial code (and requires one to write a 
full 3D code in the first place). We believe that there 
is still a need for a code able to keep the equations reg- 
ular in curvilinear coordinates while still allowing quite 
general gauge choices. 

Recently, one of us [13] presented a general procedure 
to deal with the irregularities at the origin in the case of 
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spherical symmetry. Such procedure essentially consists 
in the introduction of auxiliary variables which allow one 
to impose all the required regularity conditions on the 
metric coefficients. This method, however, cannot be 
easily extended to the case of axial symmetry without 
spoiling the hyperbolicity of the system evolution equa- 
tions. 

In this paper we follow the idea presented in [ll|, , 
and we use the general form of the tensor components 
in an axially symmetric spacetime to show that one can 
develop a generic algorithm for regularizing the evolu- 
tion equations in both axial and spherical symmetry. We 
start by writing the general form of the spatial metric 
with the corresponding symmetry. After that, we ana- 
lyze the different conditions that the geometric variables 
must satisfy at the origin or the axis of symmetry. These 
conditions arise both from parity considerations and lo- 
cal flatness. We then introduce new variables as com- 
binations of metric components whose parity properties 
guarantee that both types of conditions are satisfied at 
the same time, and evolve those variables instead of the 
original metric components. 

This paper is organized as follows. In Sec. |TT] we 
present the dynamical variables and the evolution equa- 
tions, both for the ADM system and for a strongly hyper- 
bolic formulation. In Sec. lIIII we introduce the regulariza- 
tion procedure for the particular case of spherical sym- 
metry. Later, in Sec. llVl we generalize this regularization 
procedure for the case of axi-symmetric spaces. In SecFVl 
we show some numerical examples in both spherical and 
axial symmetry. We conclude in section IVll In addition, 
in Appendix|X]we show explicitly that our equations sys- 
tem, in the spherical case, is manifestly regular. 



ated with gij, K = g'^^Kij the trace of the extrinsic 
curvature and Rij the three-dimensional Ricci tensor. 
In the above equations we have introduced the notation 
d/dt := d/dt — Cp, with the Lie derivative with re- 
spect to the shift. 

These evolution equations are subject to the Hamilto- 
nian and momentum constraints, which in vacuum take 
the form 



H := R- K,jK'^ ^0 , 



(2.3) 
(2.4) 



It is now well known that the ADM system of evolu- 
tion equations presented above is not strongly hyperbolic. 
However, the question of well-posedness is independent 
from the issue of the regularity of the evolution system at 
the axis of symmetry, and in what follows we will come 
back to the ADM system in order to show that the regu- 
larization procedure proposed here will work for arbitrary 
formulations of the evolution equations. However, we will 
also introduce below a strongly hyperbolic system. 



B. Hyperbolic evolution system 

Having a well-posed system of evolution equations is 
crucial in order to have a successful evolution code. Many 
different well posed formulations of the 3-1-1 evolution 
equations have been proposed in the literature. For sim- 
plicity, here we will follow the work of Nagy and collabo- 
rators , but will adapt it to the case of axial symmetry. 

We start by defining the new dynamical quantities 
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flat ) 



(2.5) 



II. EVOLUTION EQUATIONS 

Since we are interested in finding a regularization al- 
gorithm that is generic in the sense that it can be used 
with any formulation of the evolution equations, we will 
introduce here two different systems of evolution equa- 
tions as test cases, namely the standard ADM system 
and a strongly hyperbolic system. 



A. ADM evolution system 

We will start from the standard ADM evolution equa- 
tions in vacuum 



— gtj = ~2aK,j , 
d 

dt ^ ^ 



(2.1) 



R,,-2KuK',+KK,A, (2.2) 



where a is the lapse function, the shift vector, gij 
the spatial metric, Vi the covariant derivative associ- 



with r5^„ the Christoffel symbols associated to the metric 
gij in some curvilinear coordinate system, and rj„„ jfiat 
the Christoffel symbols for flat space in the same coordi- 
nate system. As already mentioned in fl^ , the quantities 
AJ„„ are components of a well defined tensor, while the 
rj„„ are not and in fact are not even regular in spherical 
coordinates. One must also remember that the contrac- 
tion used to construct the vector A' = g'""Am„ must be 
done with the full metric associated with the space under 
study, instead of the flat metric. 

We will now promote the A* to independent variables. 
Using the evolution equation (|2.ip we flnd the following 
evolution equation for the vector A' 



d_ 
di 



a (2K™ - 



+ 2 a if'" A 



Im 



(2.6) 



In order to study the hyperbolicity of the system, we 
must also say something about the evolution of the gauge 
variables a and /3*. For the lapse, we will choose a slicing 
condition of the Bona-Masso family |T3| of the form 



dt 



a — —a^f{a) K , 



(2.7) 
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where /(a) is a positive but otherwise arbitrary function 
of a. We will also assume that the shift vector is an 
a priory known function of spacetime x^), so that 
its derivatives can be considered as source terms for the 
hyperbolicity analysis. 

Since we want to work with a first order system of 
equations, we define the auxiliary variables 

AjTc := ^ digjk , Fi := d^lna. (2.8) 

The evolution equations for Fi and Dijk can be ob- 
tained directly from (|2.ip and (|2.7p . Up to principal part 
these evolution equations take the form 



(2.9) 
(2.10) 



where now do = dt — and where the symbol ~ indi- 
cates equal up to principal part. 

In order to obtain a hyperbolic system, we will also 
modify the evolution equations (|2.6p for the vector A; 
by adding to them a multiple of the momentum con- 
straints ([131): 

do A, ~ -a (2 ~ +2aM, 

~ -adiK . (2.11) 

For the evolution equations for the extrinsic curvature 
Kij we start by writing the Ricci tensor as 

R^J = -^5'' dAg,, + gki^^,) {A'= + r'^ |fl,t} 

+ g'"r^,rfe™, -f {A' + r'iflat} (2.12) 

+ i g'^'g"" { 2 d(igindmgj)k - d(^,gindj)gkra] , 

where F' |flat= 9^'' ^Ij I flat- In the last expression the 
symmetrization refers to the i,j indexes only. The prin- 
cipal part of Ricci tensor becomes 
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g didmgij + 5(,Aj) 



- 9(,A,) . (2.13) 

The evolution equation for Kij can then be written as 

doK,,^-adkk% , (2.14) 

where we have defined 

A*;,. :=i?% +5fjF,)-A,)] . (2.15) 

Our system of evolution equations then takes the final 
form 



doF, ~ -afd,K, 

doDijk ~ -a diKjk , 

doAi ~ -a diK , 

doK,, ~ -adkAl^ . 



(2.16) 
(2.17) 
(2.18) 
(2.19) 



Even though the A^^- are not independent quantities, it 
is very useful for the subsequent analysis to write down 
their evolution equations. Using (|2?T6l) . ([2l7| and ([2?T8l) 
we find 



doAf 



g'' diK.j 



if 



l)5%d,)K 



(2.20) 



We then have a system of 30 equations to study, corre- 
sponding to the 3 components of Fi, the 18 independent 
components of Dijk, the 6 independent components of 
Kij, and the 3 components of A^. To proceed with the 
hyperbolicity analysis we will choose a specific direction, 
say X, and ignore derivatives along the other directions. 
The idea is then to find 30 independent eigenfunctions 
that will allow us to recover the 30 original quantities, 
where by eigenfunctions here we mean linear combina- 
tions of the original quantities u = {Fi,Dijk, Kij, Ai), of 
the form Wa = J2b C-abUb, that up to principal part evolve 
as dtWa + XadxWa — 0, with Xa the corresponding eigen- 
speeds. 

Taking then into account only derivatives along the 
X direction we immediately sec that there are 14 eigen- 
functions that propagate along the time lines with speed 
— /?^, namely Fq and Dqij for q ^ x. Furthermore, 
taking / times the trace of (|2.17|) and subtracting it 
from (|2.16p . we find that the 3 functions Fi — JDim^ also 
propagate along the time lines. Finally, subtracting the 
trace of ([2T7l) from ([218]) . we find that the A""™ - A^ 
are 3 more eigenfunctions that propagate along the time 
lines. Thus, we end up with 20 eigenfunctions propagat- 
ing along the time lines with speed —f3^. 

The remaining 10 eigenfunctions are obtained by com- 
bining the evolution equation for the extrinsic curva- 
ture (|2.19p with the evolution equation for the A'^y , equa- 
tion (|2.20p . For simplicity, we assume that (3^ — Q. 
Therefore, ii q ^ x we obtain the system 



doKqi 
5oA?, 



-a d^K^, , 
-ag''''dxK, 



qi I 



(2.21) 
(2.22) 



from which is clear that we have 8 new eigenfunctions of 
the form 



with characteristic speed given by ±a - 



(2.23) 

^. Finally, tak- 
ing the trace of the extrinsic curvature and of A^^ , we find 

doK ~ -a 9, A" , (2.24) 

9oA^ ~ ^afg'^^dxK, (2.25) 

So that our final pair of cigcnfunc- 



with A" :— g 
tions are 



mn A X 

mn 



(2.26) 



VT^if ta^ , 

with characteristic speed zta V / g^^- 

In this way we see that for the evolution system where 
the vector A' has been promoted to an independent vari- 
able, and a multiple of the momentum constraint has 
been added to its evolution equation, one can obtain a 
complete set of independent eigenfunctions, showing that 
the system is indeed strongly hyperbolic. 
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III. REGULARITY IN SPHERICAL 
SYMMETRY 

A. Parity conditions 

There are in fact two types of regularity conditions 
for the metric components. One set of conditions comes 
directly from symmetry considerations. In spherical sym- 
metry we can write the metric quite generally as 



- grrdr^ + gee dfl'^ , 



(3.1) 



where a, (3^, and ggg are functions of r and t only, 
and dfi^ is the solid angle element: d^^ — d9 + sm9 dcj)^ . 
Spherical symmetry means that a reflection through the 
origin should leave the metric unchanged. By making the 
transformation r —r in the above metric we see that 
this implies that 



a{r) 

9rr{r) 

9ee{r) 



a{-r) , 

9rr{-r) , 

geei-r) , 



(3.2) 
(3.3) 
(3.4) 
(3.5) 



or in other words, a, g„ and ggg must be even functions 
of r, while (i"^ must be odd. The parity of the spatial 
metric coefficients clearly must be inherited by the cor- 
responding components of the extrinsic curvature, so that 
Krr and Kgg must also be even functions of r. 



B. Local flatness 

Parity considerations are not enough in order to have a 
regular evolution. There are extra regularity conditions 
that the geometric variables {gij , Kij ) have to satisfy at 
the origin that are a consequence of the fact that the 
manifold must be locally flat. 

Local flatness implies that close to the origin one 
should be able to write the spatial metric as 



(3.6) 



with f a radial coordinate that measures proper distance 
from the origin. If we now change the radial coordinate 
to some new coordinate r related to f through f = r{r), 
the metric will transform into 



dr 



(3.7) 



Expanding now f in Taylor around the origin we find 

' dr^ 



dr 



so that close to the origin we will have 



(3.8) 



(3.9) 



In other words, for any arbitrary radial coordinate r the 
metric at the origin must be proportional to the flat met- 
ric {i.e. it must be conformally flat). Taking this result 
together with the parity conditions derived in the last 
section we see that we can rewrite the spatial metric in 
spherical symmetry as 

df = Adr^ + r^TdVt^ , (3.10) 

where A and T are such that close to the origin 



T = To 



(3.11) 



with Aq = To functions of t only. 

The results just described where in fact already pre- 
sented in In that reference the condition that 
Ao{t) = To{t) was implemented by defining a new dy- 
namical variable that is odd at the origin: 



1 f A 



(3J2) 



and deriving an evolution equation for it. Such a regular- 
ization procedure works well in spherical symmetry, but 
its direct generalization to the case of axial symmetry has 
one very serious drawback. The problem arises because 
such an algorithm introduces terms of the form SzA/p, 
with p and z cylindrical coordinates, that change the 
characteristic structure of the evolution equations and 
can therefore spoil the hyperbolicity of a given formula- 
tion. 

Because of this, we will introduce here a different reg- 
ularization procedure that can be generalized more di- 
rectly to the case of axi-symmetry. Let us start by defin- 
ing the variables 



H := 



A + T 



J := 



A~T 
2r2 



(3.13) 



The results derived above imply that both H and J are 
regular functions that are even at the origin. The defini- 
tions of H and J can easily be inverted to give 



A-.^H + r'^J, T:=H-r^J, 
so that the spatial metric can be rewritten as 



H 



^J) dr^ (H-r^j) dn^ 



(3.14) 



(3.15) 



In order for this form of the metric to be maintained 
in time, one must ask for the extrinsic curvature to be- 
have in the same way. We will then take the extrinsic 
curvature to be 



Ka 
Ki, ^ \ t^Kt 

r^ifysin^i 



where Ka = Kh + r^Kj and Kt = Kh — r 
with {Kh,Kj) even functions at the origin 



(3.16) 



'^Kt. and 
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The A' vector in this case takes the simple form 
= (A''(t,r),0,0), where 



A" 



I Dr 



A \ A 



2 {Dree-2rJ) 
T 



(3.17) 



In this last expression we used the definition (|2.8p for the 
spatial derivatives. The parity properties of A'' follow 
directly from those of the metric, and one finds that A*" 
must be odd at the origin. 



C. Regularization algorithm 

The main idea of the regularization algo- 
rithm is simply to evolve directly the variables 
{H, J,drH,drJ, Kh, Kj, A^) imposing the appro- 
priate parity conditions on these variables, which will 
automatically guarantee that local flatness is maintained. 

The parity conditions are in fact very easy to im- 
plement numerically. The easiest way to do this is to 
stagger the origin, with a fictitious grid point located at 
r = —Ar/2. One then implements the parity conditions 
across the origin by simply copying the value of a given 
variable from r = Ar/2 to r — —Ar/2, with the appro- 
priate sign. 

The evolution equations for H, J, drH and drJ are in 
fact trivial to obtain. For example, in the case of zero 
shift, they have the form 



dtH = 


-2aKH , 


(3.18) 


dtJ = 


-2aKj , 


(3.19) 


dtDn = 


-2dr(^aKH) , 


(3.20) 


dtDj - 


-2dr[aKj'^ , 


(3.21) 



where we defined drH = Dh and drJ = Dj. The evo- 
lution equations for Ku and Kj can also be obtained 
directly from those of Ka and Kt- The resulting equa- 
tions are again trivial to derive but rather long, and we 
will write them explicitly in the Appendix |^ However, 
the evolution equation for Kh looks like 



dtKH 



aH^ 



2rA^ T2 



A" H^ - FrH -2D 



H 



(3.22) 



where Ti. stands for terms that are not divided by r. By 
simple inspection one can see that all terms in the evolu- 
tion equation for Kh are manifestly regular. The evolu- 
tion equation for Kj, on the other hand, takes the form 



dtKj = 



a H^ {h A"" - Fr^ a H^ {h drA'' - drFr 



2 r3 T2 



2 r2 T2 



J 



(3.23) 



where stands for terms that cither have no divisions 
by r, or else involve terms of the form {Dfj)"^ /r"^ , Dj/r, 
etc., which are manifestly regular. 



In the above equation, one must remember that be- 
cause of the behavior of A'' and F^., drA^ and d^F^ are 
even functions at the origin. One can now see that the 
first two terms in (|3.23p are regular by first noticing that 
they can be joined in pairs to form a single derivative, so 
that the equation becomes 

„ ^ aH^ ^ /A^\ aH^ „ / Fr 



2rA^ T2 



J 



(3.24) 



It is now easy to see that this last evolution equa- 
tion is manifestly regular, due to the fact that 
A' /r - constant + 0{r^), so that dr {A^'/r) - 0(r), and 
Fr/r ^ constant + 0{r^), so that dr{Fr/r) ~ 0{r). 

One can also see that the evolution equation for A'', 
and both the Hamiltonian and momentum constraints, 
are trivially regular. On the other hand, if one uses 
the regularization procedure of [lij , the momentum con- 
straint remains irregular. 



IV. REGULARITY IN AXIAL SYMMETRY 

A. Parity conditions 

In the case of axial symmetry, the spacetime metric 
can be written in cylindrical coordinates (p, z,(j)) as 

+ Qppdp^ -t- Qzzdz'^ -I- g^^d(j)^ 

+ 2{gpzdpdz + gp^dpd4> + gz4,dzd4)) . (4.1) 

As before, axial symmetry implies that the metric should 
remain unchanged under the transformation p —p, 
which implies 

aip) - a{-p) , (4.2 

Pp{p) - -Ppi~p) , (4.3 

(3z{p) - M-p) , (4.4 

Mp) - M-P) . (4-5 

9pp{p) = 9pp{-p) : (4-6 

gzz{p) = gzz(-p) , (4.7; 

500(p) = 94>4>{~P) > (4-8 

gpzip) = -gpz{-p) , (4.9 

9p4>{p) = -9p<p{-p) ' (4-10 
5z0(p) = 9z^{-p) ■ (4.11 

Again, the components of the extrinsic curvature inherit 
their parity properties from the corresponding metric co- 
efficients. 

B. Local flatness 

As in the spherical case, parity conditions are not 
enough. One also needs to consider the conditions aris- 
ing from the fact that space must be locally flat at the 
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axis of symmetry. We will derive those conditions here 
somewhat informally in order to have a more intuitive 
idea of where they come from. For a more formal proof 
the reader can look at where the same conditions 
are arrived at by solving the Killing equation for axial 
symmetry. 

Let us start by considering the general spatial metric 
in Cartesian coordinates 

dl^ = gxxdx^ + gyydy"^ + gzzdz^ 

+ 2gxydxdy + 2gxzdxdz + 2gyzdydz . (4.12) 

Axial symmetry implies, in particular, that the metric 
must be invariant under reflections about the x and y 
axes, and under exchange of x for y. Local flatness also 
implies that the metric must be smooth. These two re- 
quirements together imply that for fixed z we must have 



gxx ~ kp + 0{x^ +y^) - 


kp - 


I-O(p'), 


(4.13) 


gyy kp + 0{x^+y^)- 


kp - 


l-O(p'), 


(4.14) 


g,, ^ k, + 0{x^ + y^)- 




KO(p'), 


(4.15) 


gxy - 0{xy) ^ 0{p^) , 






(4.16) 


gxz ~ 0{x) ~ 0{p) , 






(4.17) 


gyz ^ 0{y) ^ 0{p) , 






(4.18) 



where kp and kz are constants. Let us now consider a 
transformation to cylindrical coordinates (p, z^cf): 



X — p cos <p , y — P sm , z — z 
Under such a transformation we have 

,-,2 J, I „ -2 



(4.19) 



gpp = gxx COS (I) + gyy sin 4> 

+2gxy sin (f) cos (j) , (4.20) 

gzz = gzz , (4.21) 
= P^ {gxx sin^ (/> + gyy cos^ (p 

—2gxy sin (j) cos (j)) , (4.22) 

gpz = gxz cos (p + gyz sin (j) , (4.23) 
9p<j> = p {gyy - gxx) sin (j) cos (j) 

+pgxy{cos'^(f>~sin^(f>) , (4.24) 

gz,p = p{~gxzsin(l) + gyzcosc/)) . (4.25) 

From the behavior of the different Cartesian metric 
components near the axis we then see that 



9pp ^ 


- kp + 0{p^) , 


(4.26) 


gzz r- 


^ fc,+0(/72). 


(4.27) 


9<t><t> " 


^ p^kp + 0{p^)) , 


(4.28) 


9pz ^ 


- 0{p) , 


(4.29) 


9p<P ^ 


- O(p') , 


(4.30) 


9z^ - 


- O(p') • 


(4.31) 



Therefore the spatial metric can be written as 

df = Adp^ + Bdz^ + p^Td(j)^ + 2(^pCdpdz 



+ p^ Ci dp dcj) + p^ C-2. dz d(pj 



(4.32) 



with {A, B,T,C,Ci,C2) all even functions of p on the 
axis. Again, let us define the new variables 



H 



A + T 



J := 



A-T 



2p2 



(4.33) 



The results derived above imply that both H and J are 
regular functions that are even in p. The definitions of 
H and J can easily be inverted to give 

A:=H + p^J, T:=H-p'^J, (4.34) 

so that the spatial metric (I4.32p can be rewritten as 

dP = [H + p^J)dp'^ + Bdz"^ + p'^{H - p^J)d(t)'^ 

+ 2{pCdpdz + p^C'idpdcj) + p^C2dzd(j)). (4.35) 

For the extrinsic curvature Kij we take the similar 
form: 



Ka, pKc p^Kc, 
K,, = \ pKc Kb p^Kc, 
p'^Kc, p'Kc, p'Kt 



(4.36) 



with Ka = Kh + p^Kj,Kt = Kh - p^ Kj. The 
extrinsic curvature components are given in such a way 
that all the functions are even, as in the metric case. 

The A' vector takes the form (A'', A^, A"^), and is a 
well defined vector. The general expression for A* can be 
obtained directly from its definition. In this way, we find 
that A'' is odd, while A^ and A"^ are even with respect 
to reflections on the axis. 



C. Regularization algorithm 

The main idea of the regularization algorithm is again 
to evolve directly {H, J, DpH, DpJ, D^H, D^J, Kh, Kj), 
instead of {A,T, Dppp, Dp^z, D^pp, D^^z, Ka, Kt), to- 
gether with the other metric and extrinsic curvature coef- 
ficients and the A' . The corresponding parity conditions 
can again be implemented numerically by staggering the 
axis with a fictitious grid point located at p = —Ap/2. 

The evolution equations for Kh and Kj can again be 
obtained directly from those of Ka and Kt- The result- 
ing equations are very long so we will not write them 
here, but they are again trivial to obtain. Consider, for 
example, the case of the hyperbolic system without rota- 
tion and, by simplicity, shift vanish. That is, equations 
(lO) . (1^ . (EH) and (I23D with Ci = C2 = 0. In this 
case the evolution equation for Kh is manifestly regular, 
but the evolution equation for Kj has terms that at first 
sight appear irregular and have the form 



dtKj = - 



f HAP 



HdpAP 




(4.37) 
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where again J stands for terms that either involve no 
divisions by or involve terms like [DpH dpFp)/ , 
DpJ / which are manifestly regular. 

Just as in the spherical case, one can see that the first 
terms in (|4.37|) are regular by noticing that they can be 
joined in pairs to form a single derivative, so that the 
evolution equation for Kj becomes 



A. Minkowski in spherical symmetry 

As a first example of the regularization method we 
evolve Minkowski spacetime with a non-trivial slicing and 
vanishing shift, using the hyperbolic system presented in 
Section IIIBl The initial data corresponds to a trivial 
slice so that 



dtK 




(4.38) 



It is now easy to see that this last evolution equation is 
regular, due to the fact that Ap / p ~ constant -t-O(p^), so 
that, dp (AP/p) ~ 0{p), and Fp/p - constant + 0{p'^), 
so that dp{Fp/p) ^ 0{p). On the other hand, by inspec- 
tion one can see that all terms in the remaining evolution 
equations are manifestly regular leaving us with a regular 
system of equations for the axial symmetric case. As final 
comment, notice that since the regularization algorithm 
is very general, one can use it in order to have a regu- 
larized system in the ADM case. One obtain the similar 
evolution equation as in the hyperbolic system. For ex- 
ample, the evolution equation without rotation and shift 
vanish for Kj looks like, 



dtKj 



aBH* ( Dp,, 
2pA^T^{AB - p-^C^Y \ "y P 




(4.39) 



One can see that the last equation is regular on the axis. 



V. EXAMPLES 

In the simulations shown below we will see how the reg- 
ularization procedure described in the previous sections 
works in practice. We will consider first an evolution of 
Minkowski spacetime with a non-trivial slicin g in order 
to compare with the algorithm presented in [lj|. We 
will perform similar simulations using both a spherically 
symmetric and an axially symmetric code. Also, in order 
to see that the regularization procedure is independent 
of the hyperbolicity of the system of evolution equations, 
we will do the axi-symmetric simulation using both the 
ADM system and the strongly hyperbolic system derived 
in Section fllBI As a second example, we will consider a 
Brill wave spacetime as a non-trivial test of the regular- 
ization procedure in axi-symmetry. 

All runs have been performed using a method of lines 
with iterative Crank-Nicholson integration in the time, 
and standard second-order centered differences in space. 



which implies. 



A = T=l , 
Ka = Kt^O, 



H = 1 , J = . 



(5.1) 
(5.2) 



(5.3) 
(5.4) 



In order to have a non-trivial evolution, we chose a non- 
trivial initial lapse profile of the form: 



i(i = 0) = 1 + r^R exp 



r-ro 



exp 



r + ro 



(5.5) 



This specific lapse profile has been chosen because it 
guarantees that a{t = 0) is regular at the origin. In 
the simulation shown below we have taken the Gaussian 
parameters to be i? = 0.001, tq = 5.0 and a = 1.0. We 
will evolve the lapse using a Bona-Masso slicing condi- 
tion but restricted to harmonic slicing, that is f{a) = 1. 
Furthermore, we have used a grid spacing of Ar = 0.1 
with the outer boundaries at r — 100, and a Courant 
factor of At/Ar = 0.5. 

Figure [T] shows the evolution of the lapse function a 
near the origin. One can clearly see that the lapse re- 
mains perfectly smooth when the Gaussian pulse goes 
through the origin. The system can in fact evolve for 
very long times and the behavior at the origin remains 
well behaved. 



B. Minkowski in axial symmetry using ADM 

In this section we present a similar evolution to the one 
of the last Section, but done now with an axi-symmetric 
code using the ADM formulation. We again consider 
initial data corresponding to a trivial slice of Minkowski, 
so that the initial metric and extrinsic curvature have the 
form 



A 
C 
Ka 
Kc 



B = T =1 , 
Ci = C2 = , 
Kb = Kt = Q. 
Kc, ^Kc.^Q 



(5.6) 
(5.7) 
(5.8) 
(5.9) 
(5.10) 
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FIG. 1: Evolution of Minkowski spacetime with a non-trivial 
slicing using a spherically symmetric code. The plots show 
the evolution of the lapse function a at different times. Notice 
how smoothly the lapse behaves as the Gaussian pulse goes 
through the origin. 




■ -8 -6 -4 -2 2 4 6 8 -8 -6 -4 -2 2 4 6 8 

FIG. 2: Evolution of Minkowski spacetime with a non-trivial 
slicing using an axi-symmetric code with the ADM formula- 
tion. The plots show the evolution of the lapse function a at 
different times along the p axis. 



which imphes, 



i? = 1 , J = 
Kh = Kj = 0, 



(5.11) 
(5.12) 



Again, we chose an initial non-trivial lapse profile which 
for simplicity is now centered at the origin: 



a{t^O) = l + i?exp 



(5.13) 



For the particular simulation presented here we have 
taken R = 0.015, <Tp = = 2.5. We will again evolve 
the lapse using harmonic slicing, /(a) = 1. For this 
simulation we have used a grid spacing of Ap = Az ~ 
0.125, with a Courant factor of At/Ap = 0.25. The 
outer boundaries are at p = 129, z = ±129. Furthermore, 
Kreiss-Oliger fourth order dissipation has been added for 
stability |18j . whereby we modify a given variable u in 
the new time-step by adding to its evolution equation 

dtUij dtUtj - {ui+2,j - 'iui-ij 

+ 6Utj - 4Mi_1j + Ui-2,j) 

+6uij - 4Mij_i -I- Mi,j_2) , (5.14) 
where the indices i,j refer to the grid points along the 



p and z directions. For this simulation the dissipation 
coefficients have been taken to be ep — = 0.04. 

Figure [2] shows the evolution of lapse function a along 
the p axis. Notice that again there is no problem at 
the axis of symmetry: The lapse evolves as a wave, goes 
through the origin, and finally returns to 1. The evo- 
lution time is only limited by the instabilities produced 
from the fact that ADM is not strongly hyperbolic. 



Minkowski in axial symmetry using a 
hyperbolic formulation 



In our next example, we consider exactly the same 
situation as in the last Section but using now a hy- 
perbolic formulation. As before, we have used a grid 
spacing of Ap — Az — 0.125 and a Courant factor of 
At/Ap — 0.25. Again, Kreiss-Oliger second order dissi- 
pation has been added for stability with dissipation co- 



efficients 



0.04. 



Figures [3l|4] and [5] show the evolution of the lapse func- 
tion a, the radial metric A, and the variable A'' along the 
p axis. We evolve the system until 50M and all variables 
remain well behaved on the axis. 
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FIG. 3: Evolution of Minkowski spacetime with a non-trivial 
slicing using an axi-symmetric code with a hyperbolic formu- 
lation. The plots show the evolution of the lapse function a 
at different times along the p axis. 



FIG. 5: Evolution of Minkowski spacetime with a non-trivial 
slicing using an axi-symmetric code with a hyperbolic formu- 
lation. The plots show the evolution of the A'' function at 
different times along the p axis. 



D. Brill waves 
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FIG. 4: Evolution of Minkowski spacetime with a non-trivial 
slicing using an axi-symmetric code with a hyperbolic formu- 
lation. The plots show the evolution of the metric function A 
at different times along the p axis. 



For our last example we have considered a non-trivial 
Brill wave spacetime, which corresponds to strong non- 
linear gravitational waves in vacuum. The construc- 
tion of such a spacetime starts by considering an axi- 
symmetric initial slice with a metric of the form 

ds^ = [e^" [dp" + dz^) + pHf] , (5.15) 

where both q and ^' are functions of (<, p, z) only. In 
order to solve for ^P, we first impose the condition of time 
symmetry, that is, Kij = 0. This condition implies that 
the momentum constraints (|2.4p are identically satisfied. 
We then choose a specific form for the function q and 
solve the Hamiltonian constraint for 5*, which for the 
metric (|5.15p becomes 

+ J (9,PP + 9,-) * = , (5.16) 

with the flat space Laplacian. The function q is quasi- 
arbitrary, and must only satisfy the following boundary 
conditions 

9lp=o = 0, (5.17) 
d]^q\p=o = for odd n , (5.18) 
q\r^oo = 0(r-2) . (5.19) 

Once a function q has been chosen, all that is left for one 
to do is to solve the elliptic equation (|5.16p numerically. 
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FIG. 6: Conformal factor \E'(p, z), along of the p and z axis 
respectively, for Brill initial data with a source function q of 
the form (|5.2Up and an amplitude of o = 2. 



Different forms of the function q have been used by 
different authors 11?, '23. Here we will consider the one 
introduced by Holz and collaborators in [20J, which has 
the form 



q 



(5.20) 



with a a constant that determines the initial amplitude 
of the wave (for small a the waves disperse to infinity, 
while for large a they can collapse to form a black hole). 
Figure [H] shows the value of \1/ along the equator z = 
and axis p = obtained by solving equation (|5.16|) 
numerically for an amplitude of a = 2.0 (small enough 
so that no black hole is formed, but large enough so that 
we are far from the linear regime). 

For the evolution of this initial data we have used a 
grid spacing of Ap = Az = 0.125 and a Courant fac- 
tor of At/ Ap — 0.2, with the outer boundaries located 
at jO = 165, z = ±165. For the lapse evolution we use a 
1+log slicing condition, which corresponds to a Bona- 
Masso (|2.7|) slicing with /(a) = 2/a. Again, Kreiss- 
Oliger second order dissipation has been added for sta- 
bility. 

Figures [3 [8] and [9] show the evolution of the extrinsic 
curvature component K^z, and the metric components 
B and T along the p axis. Notice again that for this 
simulation there is no problem at the axis of symmetry 
in the evolution of the different geometric quantities. 



VI. DISCUSSION 

We have presented a regularization procedure for the 
numerical simulation of spacetimes with either spherical 
or axial symmetry following an idea of Rinne and Stew- 
art [Tlj . This procedure enforces both the parity con- 
ditions and the conditions arising from local flatness at 
the origin and the axis of symmetry. We paid particular 
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FIG. 7: Brill wave simulation. The plots show the evolution 
of the extrinsic curvature component K^z at different times 
along the p axis. 



attention to the fact that our regularization procedure is 
independent of the system of evolution equations chosen, 
explicitly showing this in the case of the ADM formula- 
tion, as well as a strongly hyperbolic formulation similar 
to that of Nagy, Ortiz and Reula fl^ (slightly modified 
in order to have all the dynamical variables well defined 
in curvilinear coordinates). 

We have also described numerical codes that follow 
such regularization procedure both in spherical and axial 
symmetry, and presented several examples clearly show- 
ing that all dynamical variables remain regular at the ori- 
gin and axis of symmetry in each case (similar numerical 
experiments using the regularized Z4 system of 111] . [l 3| 
have also been carried out by Rinne and Stewart in p^). 
These results show that one can construct well behaved 
numerical codes in both spherical and axial symmetry 
that can allow the study of interesting astrophysical sys- 
tems with quite modest computer resources by using 
symmetry adapted coordinate systems. 

We conclude by mentioning that one can also construct 
regular codes using specialized gauge choices which can 
even allow one to reduce the number of independent com- 
ponents of the metric. Nevertheless, our interest here has 
been to find a regularization procedure that works in the 
most general case, while leaving the choice of gauge com- 
pletely arbitrary. 
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FIG. 8: Brill Wave simulation. The plots show the evolution 
of the metric function B at different times along the p axis. 



APPENDIX A: EVOLUTION EQUATION IN THE 
SPHERICAL CASE 



In this appendix we give explicitly the evolution equa- 
tions for the hyperbolic system presented in the sec- 
tion IIIBI in order to show that all equations are man- 
ifestly regular. We only consider the case of spherical 
symetry, since in axi-symmetry the final equations are 
just too long to write down here (they have been calcu- 
lated with Mathematica, and some of them are dozens of 
lines long). For the spherical case, the resulting system 
of evolution equations is, together with (|3.2ip : 



dtKn 



(A''ij2 - FrH - 2Dh) 



H^Jr^ (J + 2KaKh - 2rFrJ) + {A^DhT^ - rJ {r^H^J + r^AHJT + A^T'^)) A'' 



„5 t3 



„6 73 



J* {-2KaKh + 3 J (riv - 5)) - A^F/T 



Dj {AF - H^J - br^J^) + 2H^ ( J + 2rF^ J + Ka {Kt - r^K.,)) - iA^T^drF, 



r'^HD 



2r^HJ^ (2 J (rFr + 7) + Ka {Kt - r^Kj)) 



-DnDjr^ {3H^ - 2HJr^ + TJ^r^) + A^T^Q^ A'' - AT^drDn 



- {DjH - 36 - bDjJr'^ 



(Al) 



dtKj = TT-T^ dr[ — ]- -TT-r^ dr[ — ]+ TTTT^ i^^rH - Dh) 



2rA^T^ \ r ) 2rA'^T^ \ r ) Sr^A'^T^ ' ' 2rA^T^ 

3 Pr'^ ~ 



{-6HDj + IUDh - 2HJFr + SJH^A'') + Df^ \^-5HJ 

Dh {H^ [TDj - 2JFr) + l2HJ^r - 2HJ {D., + Jr^Fr) + WjV^ + jV^ {iDj + 2JF,.)) 
AH^KaKj + HJr^ (-SZJjr^ + 4,/^ (rF,. - 13) + 4Jr {-6Dj + KaKjt)) 

4i72 (^JKa [Kh + Kt) - SJrDj + (^2r^F^ - 9 - ^^^^ ^ - 2AT^drDj 

J^r'^ {SJKa [Kh + Kt) + l2rJD., - SD^jr^ + {rFr {rFr + 1) + 5)) 
2 {A^DjT^ + J^r {-2Jr^ + T) {H^ + AT)) A'' + 2 {2H^J^r^ - J^r^) drFr 

2 J {A + H) Jr^ + AJ^r^ + AH^T) drA' , (A2) 
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FIG. 9: Brill Wave simulation. The plots show the evolution 
of the metric function T at different times along the p axis. 



2i7V {<6JKh + AK.j) + UJ^r^Kff - lOAJ^r^Kj 
AFrT {-AK.y + Kh{A + 2Jr2)) + 2DhKaT^ + 2DjKAr'^T^ - A^TdrKn 
2H (J2 {-AKrr^ + AKjr'') + AT {FrKjr^ + drKn)) +A{A + 2 J) r^TdrKj 



(A3) 



Considering the results of Section llllj we see that the 
above equations are manifestly regular at the origin. 
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